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SUMMARY 


The purpose of this study was to obtain a theoretical formulation for the pres- 
sure on the surface of an arbitrary body moving subsonically through a coitpressible 
fluid and to implement the formulation numerically. Important applications exist in 
the areas of propeller, helicopter, and wing theory. 

Farassat has derived a solution to the Ffowcs Williams -Hawkings equation that 
describes the acoustic pressure due to an arbitrary body in motion. This solution 
has been very useful in calculating noise due to propellers and helicopters. It is a 
soltt ion to the inhomogeneous wave equation using the Green* s function approach where 
sources are distributed over the surface of the body. The pressure on the surface of 
the body is assumed known. 

As is typical of integral solutions formulated using the Green’s function, the 
solution developed by Farassat must be carefully interpreted on the surface of the 
body. The proper way to obtain the governing equation for the surface pressure is to 
integrate exactly over a small region of the body, take the limit as the "observer” 
approaches the surface, and then let the size of the region vanish. This procedure 
is well-known for a doiiblet distribution and the Kirchhoff -Helmholtz equation. As 
part of this study, this procedure was carried out for Farassat* s equation, a non- 
trivial problem since the integral is over a body in three-dimensional space and 
time. The doublet distribution and the Kirchhof f-Helraholtz equation can be shown to 
be special cases of the resulting integral equation. 

The integral equation for the surface pressure is amenable to standard nximerical 
integration techniques, but special precautions must be taken in the vicinity of the 
singularity on the surface of the body. By breaking the surface into N regions 
and using quadrature over these regions, one obtains a linear algebraic equation. 
Repeating this procedure for the observer at N distinct points gives a system of 
linear equations, whose solution yields the pressure at each of the N points. 

These pressures can then be used as inputs to Farassat* s equation to find the noise 
due to the body or to find the lift and drag. 

The above technique resembles well-known panel methods of aero<^namics (e.g. , 
Smith and Hess's method) except that coitpressibility effects are accounted for 
exactly via retarded time without recourse to the Prandtl-Glauert rule. This is 
inportant since the rule is only valid for two-dimensional or axisymmetric flows. 
Since the governing (integral) equation for surface pressure has been determined, 
the present method is a direct method, in contrast to most panel methods where 
sources and doublets are distributed over the ho6y and their strengths determined 
numerically. Another advantage of this technique over more conventional panel 
methods is that it uses the pressure as the dependent variable. These quantities are 
continuous across the wake, unlike the velocity potential. Tlierefore, there is no 
need to approximate the wake’s location or assume a distribution of sources over it. 

A coitputer program has been written to inplement the above -described theoretical 
formulation for arbitrary bodies in rectilinear and/or angular motion. Numerical 
results have been obtained for the surface pressure distribution on ellipsoids, 
wings, and helicopter rotors. The agreement between these results and those from 
other theoretical techniques and experiments has been good even for lifting bodies. 
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INTRODUCTION 


In the last few decades, renewed interest in rotating blade technology, specifi- 
cally in propellers and helicopters, has reversed a decline in propeller development 
that began in the 1940 's with the advent of gas turbines- Since turbojets, and later 
turbofans, offered significant improvements in reliability, weight, and speed propel- 
lers were relegated to only a few specific applications. The low price of fuel 
allowed the pursuit of higher and higher speeds. However, the rapid increase in fuel 
costs in the early 1970’ s changed the equation for direct operating costs, and there- 
fore the effectiveness of all forms of transportation had to be reevaluated. It was 
soon realized that propellers could be more efficient than turbofans, even at rela- 
tively high speeds (up to M = 0.8 cruise (ref. 1)). 

In the future, highly swept multibladed propellers (propfans) will iirprove pro- 
pulsive efficiency by 15 to 22 percent over advanced turbofans (at M = 0.7 to 0.8), 
even after gear and turbine losses. At lower speeds, the propulsive efficiency of 
propellers is already nearly 90 percent. The propeller definitely has a place as a 
propulsion device of the future. 

Helicopters represent a conpletely different application of rotating blades. 

They have become more and more popular, but for reasons other than efficiency. Their 
main advantage is their ability to land in small areas, which often eliminates the 
need for ground transportation from the landing pad to the destination. This ground 
travel time can be significant, since most airports are far from city centers. 
Including travel from city to airport, an airplane trip from Ijondon to Paris requires 
160 minutes; the same trip by helicopter requires 70 minutes (ref. 2) . 

Helicopters satisfy many requirements of communities, businesses, and the mili- 
tary- They are being used to deliver patients to hospitals when time is critically 
important, and they assist in law enforcement, forestry, and traffic monitoring. 

Even relatively small corporations can afford helicopters to rush employees from 
place to place. However, the largest demand for helicopters comes from the military 
for rescue, reconnaissance, and antitank roles. In addition, they are performing 
many tasks formerly done by trucks, including hauling airplanes, trucks, and tanks. 

Of primary importance is their role in rescuing; during the Korean War, helicopters 
transported over 23 000 casualties, over half of which would have otherwise died 
(ref. 2). 


With renewed interest in rotating blade technology, there is a concomitant 
interest in the acoustic and aerodynamic characteristics of rotating blades. Despite 
many advances in aeroacoustics, relatively little has been accoitplished in developing 
efficient aerodynaimic methods for propellers and helicopters. The present work 
presents a conpressible aerocfynamic method that is especially useful for rotating 
blades. Specifically, an integral equation governing the pressure on the surface of 
an arbitrary boc^ is developed, and a finite element numerical procedure (panel 
method) is used to solve it. The theory is linear and inviscid, and the dependent 
variable is the pressure, which is governed by the wave equation. 



A. Motivation 


The present study was motivated by a need for effective noise prediction methods 
for moving bodies, especially rotating blades. There are several formulations for 
the noise due to moving bodies (ref. 3). The governing partial differential equa- 
tions and boundary conditions are usually reduced to integral formulations using the 
Green's function. These formulations describe the pressure in terms of integrals 
over the surface of the body, and the surface pressure is contained in the integrands 
of these equations. Therefore, detailed surface pressure data are required to calcu- 
late the noise. 

Currently the required surface pressure data can be obtained in several differ- 
ent ways. In some cases, blade element theory is used with standard aerodynamic cor- 
rections for conpressibility. In others, e3q>erimental measurements are used. What- 
ever the technique, surface pressure is usually obtained via standard aerodynamic 
methods that were developed before the advent of high-speed conputers. 

The present work shows that the surface pressure can be calculated from the same 
integral formulations that govern the acoustics. However, in this case, instead of 
having an integral representation, one obtains an integral equation, specifically a 
singular, inhomogeneous Fredholm integral equation of the second kind. 

The key word above is singular. Without this complication, there would be few 
difficulties. Singular integral equations are ambiguous unless interpreted properly. 
This interpretation forms much of this report and is discussed in detail in sec- 
tion II. In fact, the main contribution of the present work is to derive the govern- 
ing integral equation for surface pressure and to solve it for severai different 
bodies using a panel method. 

By using the acoustic formulation to determine surface pressures and thus 
eliminating the need for ad hoc aerodynamic methods or ejqjeriment, the entire noise 
prediction methodology for moving bodies becomes autonomous. Also, as progress is 
made in modeling more and more aspects of the actual flow, they can be included both 
in surface pressure calculations and in noise calculations simultaneously. In addi- 
tion, by using the same theory to predict aeroacoustics and aerodynamics, one may be 
able to develop a deeper understanding of the two phenomena. 

The acoustic formulations presented herein were originally suggested by 
Lighthill based on the acoustic analogy (refs. 4 and 5) . Ffowcs Williams and 
Hawkings (ref. 6) e3q>anded this concept to include the effects of solid surfaces by 
writing the governing equations of fluid mechanics in the form of an inhomogeneous 
wave equation (in terms of pressure or density). The inhomogeneities represent 
boundary conditions, nonlinearities, and viscous effects. 

Using the free-space Green's function, one can derive an integro-dif f erential 
equation from the acoustic analogy. This equation does not represent a solution 
per se, but for numerical purposes it makes the problem more tractable. Farassat 
(refs. 7 and 8) has reduced this differential integral formulation to strictly an 
integral formulation, which is even more amenable to numerical treatment because the 
need for numerical differentiations is eliminated. Farassat' s equation forms the 
starting point of the present analysis. 

Using the acoustic analogy approach implies that the pressure (or density) is 
the dependent variable. In linearized theory the pressure is proportional to the 
acceleration potential, as originally discussed by Prandtl in 1936 (ref. 9). Since 
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then others have used the pressure, or acceleration, method, notably Knssner 
(ref- 10) and Kbndo (ref. 11) and, more recently. Van Holten (ref. 12) and Dat 
(ref. 13). Since the present effort was motivated by the desire for efficient noise 
prediction methods, working in the pressure formulation is natural. There was no 
interest in the details of the wake, the downwash, or the trailing vortices except 
in how they affect the pressure and normal fluid velocity on the surface- A major 
advantage in using this formulation as opposed to using the velocity potential is 
that there is no need to integrate over the wake. Since the pressure, unlike the 
velocity potential, is continuous across the wake, there are no integrals over the 
surface of the wake. Of course, the effects of the wake and trailing vortices must 
be accounted for. This is discussed in section IV. 

Linearized acoustics and linearized conpressible aerodynamics are one and the 
same. They both derive from a small perturbation of the continuity and Euler equa- 
tions, which can be combined to give the wave equation. The link between acoustics 
and aerodynamics is exploited in few works as much as it is here, however. Gener- 
ally, acoustics is concerned with the signal after it has radiated from the body. 
Aerodynamics is concerned with determining the forces on the body. In the current 
work, both effects are shown to originate from the same phenomenon. Just as one can 
find the pressure on a vibrating piston using the Kirchhoff -Helmholtz (ref. 14) 
integral equation, the pressure on the surface of any (thin or slender) moving body 
can be determined using Farassat's equation. 

As alreac^ mentioned, linearized coirpressible aerodynamics is governed by the 
wave equation; for two-dimensional or axisymmetric bodies in steady, rectilinear 
motion, the problem is simplified considerably. With the Prandtl-Glauert transforma- 
tion, the wave equation can be transformed into the Laplace equation, so that one can 
use the well-known and powerful methods of potential theory. For finite bodies 
undergoing very conplicated motions such as spinning, vibrating, and translating, the 
problem is not so simple. One has no alternative but to solve the wave equation. 

In the past, coirpressibility corrections developed for two-dimensional bodies in 
rectilinear motion were used when solving for the aerodynamics of rotating blades. 

The blade was divided into sections along the span and each section treated as though 
it were in rectilinear motion with Prandtl-Glauert or Karman-Tsien corrections for 
compressibility. Although there are more appropriate techniques such as the 
Goldstein-Lock method (refs. 15 and 16), these do not apply to helicopter blades in 
forward flight and they also require corrections for coitpressibility . The method 
presented herein represents an actual solution to the governing equation, the wave 
equation. Coirpressibility , three-dimensionality, and complicated motions are treated 
together in a unified fashion. 

Using the wave equation directly means that the problem is four-dimensional, in 
space and time. In the past this was a serious obstacle, but today with high-speed 
conputers it is not. In fact, once one becomes accustomed to the notion of four 
dimensions, the physics becomes much more understandable. 

Coirpressibility manifests itself via a finite propagation speed of disturbances. 
In this work, conpressibility effects are accounted for by considering this finite 
propagation speed. The time of propagation, or the distance the signal actually 
travels, is calculated exactly. In rectilinear motion this calculation is completely 
equivalent to the Prandtl-Glauert transformation, which "stretches” the bo<^ to 
account for the actual distance the signal travels. Because the body and the signal 
are moving, the signal must travel farther from one point on the body to another. 

In the current work by using four dimensions and the wave equation, the effect is 
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accounted for exactly for arbitrary motions. Conpressibility is discussed later 
in terms of retarded time, which has also been used by Kassner (ref. 10), Rondo 
(ref. 11), Van Holten (ref. 12), Dat (ref. 13), and Morino (ref. 17). 

After the integral equation governing the surface pressure is derived, the 
numerical technique used to solve it is discussed. The conputational method can be 
classified as a finite element technique (ref. 18), but it is more accurately called 
a panel method. Throughout this work the words element and panel are used inter- 
changeably. In a panel method, the surface of the body is approximated by a finite 
number of quadrilaterals and the pressure is assumed to follow some given behavior 
over each element; in this case it is assumed to be constant over each element. 

The numerical approach used here can also be classified as a boundary integral 
equation (BIE) method (ref. 19) which has become very popular in the last few 
decades. This method has been used in several different fields including fracture 
mechanics (ref. 20), potential theory (refs. 21 and 22), structxires (ref. 23), and 
acoustics (refs. 24 and 25) . The advantage of these methods is that the problem is 
reduced from one in space to one over a hyper surf ace. 

Probably the most common exanples of BIE or panel methods are the aeroc^namic 
codes of Hess and Smith (refs. 26 and 27). Their codes were originally designed to 
solve for the velocity potential on nonlifting bodies in rectilinear motion (ref. 28) 
and later were esqpanded to include the effects due to lift (ref. 29). The governing 
integral equation is that given by Lamb (ref. 30) as the Green's function formulation 
of the Laplace equation. This is essentially a distribution over the surface of the 
hody and the wake of sources and doublets whose strengths are adjusted to satisfy the 
appropriate boundary conditions. Numerous other very effective panel methods exist, 
notably MCAIR (ref. 31), MBB (ref. 32), PANAIR (ref. 33), and NLR (ref. 34), which 
were developed by McDonnell Douglas, Messershmitt-Bolkow-Blohm, Boeing, and the 
Netherland's National Aerospace Laboratory, respectively. 

Panel methods have been very effective for aircraft configurations, but they are 
not suitable for rotating blades because they use two-dimensional compressibility 
corrections (Prandtl-Glauert, Gothert, or Karman-Tsien rules (refs. 35 and 36)) and 
they assume rectilinear motion. Furthermore, these methods are all too conplicated 
to justify using them to calculate the inputs (surface pressures) to an acoustics 
prediction program. They would be much too costly to adapt and to run. The effi- 
ciency of the present formulation supports this claim. 

Another panel method, developed by Morino (refs. 17 and 37) , does account for 
coitpressibility in terms of retarded time. This method is an inviscid velocity- 
potential formulation in contrast to the present pressure formulation. As discussed 
earlier, one must integrate over the wake when the velocity potential is the depen- 
dent variable. Furthermore, Morino' s method results in a convected wave equation 
because the frame of reference is not fixed to the undisturbed fluid. l^parently 
he foresaw mainly aircraft-type applications for this work and thus the emphasized 
rectilinear motion. The integral equation thus acquires a very complicated form 
that would seem to limit its utility in clarifying the physics of the problem (see 
eq. 1.5, ref. 17). In addition, although both the present method and MDrino's method 
(SOUSSA-P 1.1 (ref. 38)) have been programmed for the case of linearized theory and 
both could be expanded to include nonlinear effects, the present method could also be 
ejqjanded to include viscous effects (since viscosity is alrea<^ included in the 
Ffowcs Williams-Hawkings (FW-H) equation). 
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Another important point regarding the method of reference 17 concerns the treat- 
ment of the singularity in the integral equation. Tliis treatment is also referred 
to as interpreting the equation, regularizing the equation, or putting the observer 
on the surface and is performed in section II. C of the present work. Although Morino 
criticizes Widnall's method because it requires a Cauchy principal value interpreta- 
tion, Morino does not state that his method could also be classified as such. In 
fact, so could all surface singularity methods, including the present one. This 
leads to the questions of how one treats the singularity and the meaning of the 
Cauchy principal value. In this report the singularity is treated analytically and 
its relation to the Cauchy principal value is demonstrated. Morino claims to have 
done this also (using a velocity-potential formulation) , but the results seem to 
differ from those presented here. All of this is discussed in section II. C. 

The remainder of section I is devoted to presenting the FW-H equation and one of 
Farassat's solutions to it. For completeness, the linearized version of the FW-H 
equation is derived. These sections also illustrate the importance of generalized 
function theory (ref. 39). 


B. Ffowcs Williams -Hawkings (FW-H) Equation 
The Ffowcs Williams -Hawkings (FW-H) equation is 




d [T^^ H(f)] 

5x. bx . 

1 3 


( 1 ) 


where 




signifies the wave operator. 


□ ^p = — - 

5t^ 


V^P 


This equation was first derived in reference 6. Farassat also derived it in refer- 
ence 8 using a method called embedding. It represents a combination of the mass 
continuity and conservation of momentum equations, plus the boundary conditions. 

The first term on the right hand side of equation (1) behaves like a monopole. 
It represents the normal velocity boundary condition. One can envisage this as a 
distribution of mass sources. Often called the "thickness” term, it is shown sub- 
sequently to be the equivalent to the thickness terms used in linearized 
aerodynamics . 

The second term behaves like a dipole distribution and is due to the viscous 
stresses and thermodynamic pressures acting on the surface. 

The last term on the right hand side is the quadrupole term. It contains the 
nonlinear effects, such as turbulence. In Lighthill's work (ref. 5), this was the 
only term on the ri^t hand side, but it is usually considered unimportant for noise 
due to moving bodies. However, preliminary studies recently conducted indicate that 
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these terras can be significant in the transonic and high supersonic speed regimes for 
thin wedge-shaped airfoils (refs. 40 and 41). 

The FW-H equation is easily derived from the governing equations of fluid 
mechanics using the embedding procedure (ref. 8). Ihis procedure converts the 
problem into one in unbounded space. The boundairy conditions become source terms, 
so that the free-space Green's function can be used. A simple example is given in 
reference 8. Converting boundary conditions to sources is not uncommon and is dis- 
cussed in reference 42 (pp. 791-792). 

Throughout this study only the linearized version of the FW-H equation is used; 
thus, for completeness this version is derived. The small perturbation forms of the 
governing equations of fluid mechanics are 


JL ^ 

2 5t 




o 5x. 
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(2) 


5u^ 
^o 5t 


+ 



(3) 


2 

p = c p 


where all terms higher than first order are neglected and the summation convention is 
iitplied. Tliese are also referred to as the linear acoustic equations. Iheir range 
of validity is discussed in reference 14. 

If the derivatives are interpreted as generalized derivatives (refs. 39 and 43) 
and the field variables are defined as generalized functions, then the possibility of 
finite jumps in these variables (e.g. , across surfaces and shocks) must also be 
included in these equations. If one considers only subsonic flow, the only discon- 
tinuity surface is the ho^ and the wake. It is easily shown that the governing 
equation for the pressure has no contribution from the wake. If one is interested 
in the velocity potential, the wake must be considered. In this report the wake is 
not even included because the result would be unchanged if it were. 

Interior to the bocfy the acoustic quantities are assumed equal to zero. This 
is arbitrary, however, and for other problems it may be advantageous to assume some 
other value. T^ie procedure is carried out formally by writing all the flow quanti- 
ties multiplied by a Heaviside function H(f) where f(x,t) describes the body sur- 
face and f < 0 inside the ho6y , As examples, 9p/5t and dp/dx^^ become 
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where the bars over derivatives denote generalized differentiation* Also, the 
brackets indicate the jump in that quantity across the surface. Note that inside the 
t>o(^, pressure, density, and fluid velocity perturbations all vanish. The operators 
of equations (2) and (3) are applied to the generalized quantities as follows: 
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The terms in parentheses vanish because of equations (2) and (3). Using the relation 
(ref. 30) 
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where is the velocity of a point on the body, and noting that the normal to the 

body is 

_ 1 ^ 

"i |Vf I 


gives 


dt 


-V 

n 


|Vf| 


Since there is no flow through the body. 
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Using these relations, equations (4) and (5) become (after eliminating second order 
terms ) 
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2 5t 


5u. 

X 

^o dx. 

X 


PoV^lVfl 6(f) 


( 6 ) 


7 



(7) 





5(f) 


Taking the generalized derivative of equations (6) and (7) with respect to and 

t, respectively, and subtracting yields 


□ ^P = I 6(f)] - I 5(f)] 

i 


which is the linearized FW-H equation. Terms involving products of small perturba- 
tion quantities have been neglected, along with viscous effects. It is well-known 
that near stagnation points the assumption of small perturbations is not valid, but 
this affects only a small region around the stagnation point (ref. 44, p. 209). 

An important point to make at this time concerns the pressure formulation versus 
the velocity-potential formulation. The present work uses the pressure as the 
dependent variable and the flow is linearized from the beginning. In the past it has 
been more common to use the velocity potential, especially for very low speeds, 
because at speeds much less than the speed of sound one can assume that the flow is 
incoitpressible. In this case the continuity equation becomes 


V . = 0 


or 


V 


2 


4 > = 


0 


where 


V(|) = 


Uo. 


and (() is the velocity potential and is the net fluid velocity. Thus the 

velocity potential is governed by the well-known Laplace equation. More importantly, 
however, there is no need to assume small perturbations. Thus, for very low speeds 
and inviscid flow, the problem is relatively easy to solve with all the nonlinearity 
contained in the Bernoulli equation. 

On the other hand, if one begins with the small perturbation theory of equa- 
tions (2) and (3), the governing equation is the wave equation (on the pressure, 
velocity, and velocity potential). But at very low subsonic speeds, using the wave 
equation on the pressure is not as accurate as using the Laplacian on the velocity 
potential (especially near stagnation points) , because the pressure satisfies the 
Laplacian only for small perturbations. The momentum equation is nonlinear, so that 
the Bernoulli equation is nonlinear. 
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The relation between the above theories is somewhat conplicated because it 
involves two different asyirptotic expansion procedures. In the small perturbation 
method one assumes that 


= p + ep* + . . . 


+ ep' + . . . 


U^ = £u' + ... 

where p^ and p^ are the net density and pressure, p^ is the pressure of the 
undisturbed medium, the primes indicate perturbation quantities, and e is a measure 
of the body thickness. This method becomes exact as the thickness approaches zero. 

It is often useful away from stagnation points even for relatively thick bodies. In 
the incorrpressible problem one assumes that e is a measure of the Mach number. 

This is exact as M > 0 for any body and is known to be useful for small but finite 
Mach numbers (ref. 35). 

However, since the main concern here is in rotating blades, which are thin and 
usually operate at high speeds, the linearized small perturbation theory is appro- 
priate. This formulation will be accurate as long as the disturbances remain small, 
even at relatively high subsonic Mach numbers. In fact, this formulation gives use- 
ful results for supersonic motions also. However, one must then account for multiple 
emission times. The transonic regime is, of course, inherently nonlinear. 


C. Farassat*s Solution to FW-H Equation 

Farassat has derived several different integral representations of the FW-H 
equation. Each one is particularly well suited to a different application or numer- 
ical solution technique (ref. 45). The one that is most appropriate for the present 
work is equation (9) of reference 8; 


-V 

4n p(x,t) 



f=0 


p cv + p cos 0' 
^o n ^ 

r| 1 - M I 


'ret 



dS 


( 8 ) 


ret 


where 


cos 9 = n % r 


r 


r 

r 
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Mr = 


V 


r = |x - y(T) I 


T = t 


£ 

c 


The subscript ret signifies that the expression is to be evaluated at retarded 
time i;« This accounts for cortpressibility of the flow where the "source" at y{x) 
(in motion) emits a signal that arrives at the stationary "observer" at 1c a short 
time after it is emitted. This is unlike incompressible flow where signals travel 
with infinite speed. Therefore the integration over the body surface, f(x,t) = 0, 
is not carried out at a single emission time t because different points on the sur- 
face have different emission times. Of course, one is only interested in the signals 
that arrive at x at the same time. Note that retarded time is a more common con- 
cept in electromagnetic theory than in acoustics or aerodynamics (ref. 46) . 

It is sometimes useful to consider the retarded time process in reverse order in 
terras of a collapsing sphere (ref. 8), as Nystrom and Farassat (ref. 45) did in the 
calculation of the noise due to high-speed propellers. 

It should be remembered that equation (8) is not a solution per se, but more 
accurately a representation. Since the integrals contain the unknown, the pressure 
on the surface of the bo<^, the problem has not actually been "solved." However, if 
one knows the surface pressure on a given body, equation (8) predicts the noise due 
to that body. 

Equation (8) represents the starting point for this report. Although it has 
been effective in predicting noise from bodies in motion, it has not been used to 
predict surface pressures on arbitrary bodies, mainly because the integrals become 
singular when the observer is on the surface. A mathematical limiting procedure is 
required, the details of which are described in the next section. 


II. ANALYSIS 

A. Singular Integrals and Boundary Solutions 

Singular integrals are very common in mathematical physics because of the 
frequent use of the Green's function technique, which can be thought of as a distri- 
bution of singularities or delta functions. Whenever one desires surface information 
from such a method, singular integrals may arise. These singularities must be inter- 
preted properly to obtain meaningful results, because truly singular integral equa- 
tions are ambiguous. Usually, however, these integrals are special cases of regular 
integrals, and their proper interpretation can be inferred from the physics. 

Singular integrals are especially common in aerodynamic theory because of the 
emphasis on surface data. They appear in equations for the velocity potential, the 
acceleration potential, and the downwash integral, to name a few. They are the basis 
of the modern panel methods as well as lifting surface and lifting line theories. 
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Thus the real difficulty in using the acoustic formulation to predict surface 
pressures is that the integral equation becomes singular if naively interpreted. 

That is, if x is on f = 0, then at some point, x = y, which means that r = 0. 

At this point the integrand becomes infinitely large. 

The level of difficulty that this fact presents cannot be minimized. Despite 
tremendous advances in numerical methods. Baker (ref. 47) states: 

The inpression I intend to convey is that the treatment of 
singular integral equations is not conpletely understood at the 
present.. . .There appear to exist very few effective methods for 
solving such equations. 

The proper way to interpret singular integrals is through a limiting process 
(refs. 42 and 48 to 50). In the acoustic formulations, one must assxmie that the 
observer approaches the surface in the limit from the proper side of the surface 
(inside or outside depending on the problem being solved). However, this limit must 
be taken after integrating over the surface, which presents a problem for coitplicated 
integrals such as are present here. 

As an example of how a singular integral is interpreted, a one-dimensional equa- 
tion called the downwash integral (ref. 50) is developed. The approach to interpret- 
ing this simple equation is essentially the same technique that is used to treat 
Farassat's equation. The downwash integral is defined as 


F(x,y) 


f(S) (X - 1) dz 
•'a (x - C) ^ 


(9) 


where F(x,y) represents the velocity in a two-dimensional flow due to a distribu- 
tion of sources and sinks along a < ^ < b of strength f(^). When y = 0, this 
integral reduces to the familiar integral 


f(g) dc 

Ja ^ ^ 


which is commonly defined in terms of a Cauchy principal value technique; that is, a 
small symmetric region about ^ = x is removed. Quite often this form, with y = 0, 
is the starting point for an analysis. Ordinarily the Cauchy principal value inte- 
gral, like other singular integrals in higher dimensions, derives from an integral 
that is not singular, such as the one shown for F(x,y). The surface representation 
is valid only in the limit, even for the Cauchy principal value, as is now shown. 

Strictly speaking, equation (9) is undefined for y = 0 when a < x < b. The 
proper way to obtain the value of this integral for y = 0 is by taking the limit as 
y 0, but the limiting process must be performed after the integration, that is, 

F(x,0) =lim l(i ) (X - ^) 

y^O •'a (x - 5) + y 
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Mangier (ref. 50) showed how to derive the Cauchy principal value form, by breaking 
up the region of integration. Ihat is, write 


F(x,0) 


lim 

y -^0 




f(g) (X - C) dg 
2 2 

(X - 5)^ + y^ J 


where and £2 arbitrarily small and of the same order. Now for the regions 

that do not include ^ = x, the integral is well behaved and one can set y = 0. 
Within the e-region, ejqpanding f(^) in a Taylor series about ^ = x gives 


F(x, 0) 




X - c 


lim 

y^O 


f(x) 

2 



+ 0(e) 


Taking the limit yields 


F(x,0) 


•^a -/x+e^ 


f(^) . 

X - 5 


f(x) 



0(e) 


The well-known Cauchy principal value is for an e-region symmetric about ^ = x; thus 
for e^ = £2 = £/ 


F(x,0) 




f(^) 

X - Z 


+ 


0(e) 


which is the Cauchy principal value. 

This integral is correctly termed a serai convergent integral (ref. 48); that is, 
for a differently shaped e-region, a different form of the equation is obtained. For 
a fully convergent integral, a term like ln(e 2 /e^) would not be present and thus 
the region that was given special treatment would not matter. This point is very 
important in the subsequent interpretation of Farassat's equation. To fully appre- 
ciate and use the principal value concept, one must realize that the region around 
the singularities cannot simply be deleted. Different e-regions require different 
representations. The numerical values of the final results for two different 
e-regions are the same though, because any portion the solution left out of the 
integral shows up in "extra" terms. 

In the following section an equivalent procedure is applied to Farassat's equa- 
tion to obtain an equation valid on the surface of the body. It too is semiconver- 
gent and the form of the extra term that comes from the e -region depends on the size 
and shape of the e-region, just as for the downwash integral of the Cauchy principal 
value. 
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B. Integral Equation Without Derivative 


Farassat's equation (8) is an integro-diff erential equation. Taking the deriva- 
tive inside the integral produces an additional singular integrand, that is, one with 
1/r^ dependence. Since the main interest here is in the form of the equation on the 
surface, it is important to bring the derivative under the integral to illuminate the 
singular term and allow it to be regularized. An integral equation with the deriva- 
tive inside has been used by Woan and Gregorek (ref. 51) for noise prediction and is 
presented in a more general form by Farassat (refs. 3 and 7). 

The relations necessary to eliminate the derivative are, from Farassat, 


^ = 1 ^ 
5t 1 - M 

r 



m 

5t 


• c 2 2, 

Mr + ^(Mr - M ) 


These give 


a 

1 ^ 

rM + c(M - M' 
r r 

5t 

r(1 - M ) 
L r 

”2 2 
r (1 - M ) 


In addition, one other relation is required. Farassat uses the general stress 
term lir^ and does not simplify it to the inviscid form pnj^r^ = p cos 0. This 
simplification gives 


gives 


cos 0 = 


v^ 


= r • ((A) X n) - — 


„ v^ cos 0 

n + ^ 

r r 


Defining 


Q = r 


(o) X n) 


V V cos 0 

^ on ^ I ^ 

cos 0 = Q - — + 

r r r 
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The above relation allows one to write equation (1) without the derivative. It 
becomes 


4n p( 


ip 


• 2 ^ 
(p CV + P cos 0)rrM + c(M - M )1 


2 3 

r (1 - M 

r 


dS 


ret 


/ 

f=0 


p(r^ + V cos 0 
r r 




2 2 
r (1 - M ) 
r 


dS 


ret 


r p cos 0 

f-O - ”rl 


dS 


ret 


A slightly more general form, which includes unsteady pressures and velocities, is 
derived in reference 7. 


Regrouping these terras into those that do and do not become singular on the 
surface gives 


4n p(: 


x,t) = J K (x,t;y,T) dS + f K (x,t;y, 
f=0 ^ f=0 ® 


X ) dS 


( 10 ) 


where 


^R = 


2 • • 

pcMM +pQ(1-M)+pM cos 
^o nr r r ^ r 

cr ( 1 - ^ 


/ret 


= 


p C^M (M - M^) + p[(1 - M^) cos 0 " (1 - M )M ] 
^o n r r n 

2 3 

r (1 - M^) 


ret 


For incompressible flow, is clearly a regular integrand because when x is on 

the surface f = 0, 


- dS = 0(1) 
r 


Similarly, in incoitpressible flow. Kg represents the singular portion of the 
integrand because when i is on f = 0, 



and r = 0 (x = y) at some point. 
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C. Integral Equation Valid on Body Surface 


In section II. A, a method was described for interpreting a simple semiconvergent 
integral. The final result was the well-known Cauchy principal value. In this sec- 
tion the same procedure is used to interpret Farassat's equation when the observa- 
tion point X is on the surface f(x,t) = 0. The equations are lengthy because of 
the coitplicated four-dimensional integrands, but the procedure is corrpletely analo- 
gous to that already described. As for the downwash integral, Farassat’s equation is 
undefined for r - 0 unless it is interpreted properly. 

In the case of the downwash integral, the region of integration was divided into 
one that contained the singularity and one that did not. The integral over the 
region without the singularity was numerically straightforward to compute. This was 
not true of the region that contained the singularity, subsequently called a ”hole." 
Therefore, over this hole the integral was treated analytically in such a way that 
the error could be made arbitrarily small by letting the size of the hole become 
arbitrarily small. That is, the error was shown to be of the order of the size of 
the hole. More importantly however, the form of the equation (the extra term) was 
shown to depend on the shape of the hole. Thus, the integral was semiconvergent. If 
an asymmetric hole around the singularity were chosen, the Cauchy principal value 
technique would not apply. For an asymmetric hole one cannot simply neglect the 
value of the integral over the hole because it has a finite value. This extra term 
becomes zero when the hole is symmetric. 

The first step is to divide the region of integration into two parts, one that 
includes the singularity and one that does not. For the downwash integral this was 
done by breaking the ^-axis into three parts. Farassat* s equation is an integral 
over a surface, so the surface is divided into two parts: one part over the original 

body surface with a small hole around the singularity removed and the other part over 
the surface of the hole itself. See figure 1. The combination of these two regions 
is the original surface, f(^,t) = 0. 


n 



Figure 1.- e-region of arbitrary body. 
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Thus, Farassat's equation is written 


lit p(x,t) = f Kp(x,t;y,T) dS + ^ K_(x,t;y,T) 


+ f K-(x,t;y, 


t) dS 


where / indicates that a specific hole has been removed and is the surface of 

the hole. Note that it is not necessary to break up the region of integration in the 
regular integral Kj^. 

The first two integrals present no difficulties numerically since they are both 
convergent. By definition, neither of them contains a singularity. The third inte- 
gral is difficult to calculate numerically. The proper way to obtain its form, for 
an observer located at x^, is via a mathematical limiting process where the observer 
is located at x at a distance of 6 above the surface. Then, after integrating, 
the limit as 6 0 is performed (see fig. 1). This limit is analogous to the limit 

as y 0 in the downwash integral. Note that for convenience 1c is taken to be 
along the normal to the surface at 1c^. Mathematically one can write this procedure 
as 

■> 

p(x ,t) = lim p(x,t) 

° 6-»-0 

Hence the integral equation becomes 


. r » ^ c , , 

lx p(x^,t) = J K^(x^,t;y,x) dS + j Kg(x^,t;y,' 
+ lim I K (x,t;y,T) dS 

o * — J * 


rT) dS 


6^0 f =0 
e 


The first two integrals can be obtained by simply replacing x by x^ (setting 
6=0) since they are convergent, but this is not possible in the third integral. 

One has no alternative but to perform that integration first and then take the limit; 
otherwise the integral is divergent. However, this integral is too complicated to 
allow analytic evaluation for most bodies of interest, especially since the integrand 
Kg contains the unknown pressure p. One would have to solve the entire integral 
equation and then take the limit as 6 0, to obtain an analytic expression for the 

surface pressure. Note that in equation (9), the downwash integral example, the pro- 
cedure leading to the final result was as follows; e>q>and f(^) in a Taylor series 
around ^ = x and obtain an approximate analytic expression for the integral with an 
estimate of the error. 

An analogous approach can be used on Farassat's equation. By assuming that the 
size of the hole, or e-region, is small, one can approximate the region as planar. 

In addition, the pressure can be expanded in a Taylor series about the point x^; 
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that is, the pressure is assumed to be constant over the e-region. Using these 
approximations, one can then calculate the third integral in equation (11) analyti- 
cally and obtain an estimate of the error involved. 

The remainder of this section is devoted to obtaining the appropriate analytic 
e3q>ression for the third integral in equation (11). For reasons that will be obvious 
later, a square is chosen for the shape of the e-region, but the final result is 
applicable to differently shaped e-regions with a symmetry to be described later. 

To perform this integration over a small, square, planar panel, one must first 
write the integrand Kg in terms of surface coordinates. This is nontrivial because 
the integrand involves four dimensions (space and time) . Figure 2 is an enlarged 
view of the e-region shown in figure 1. It shows the panel at two different times- 


/- ^ r 



Figure 2.- Geometry near singularity in retarded time. 


The time t is the reception time and t is the emission time. Ihus, a signal 
emitted from the surface point y(T) at time x is received at the (stationary) 
observer point x at time t. The 6 on this figure is the Scime 6 that was used 
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on figure 1. The upper plane is the position of the panel at the time t, the time 
when the observer receives the signal from y(T)» The lower plane is the position of 
the panel when the signal was emitted. The angle between the surface normal n and 
the radiation direction ^ is 0. For convenience, the coordinate system (fixed to 
the panel) is aligned so that the panel velocity vector is in the y^y 2 “plane. The 
angle between the velocity and the y^-axis is defined as <(). 

As mentioned earlier, the observation point x is always stationary. There- 
fore, when one speaks of the observation point being on the surface, one must specify 
at what time this occurs. In figure 2, notice that x > when 6 0. This 

surface point x^ is of course moving, but when the time equals t, it is a distance 
6 away from x in the direction of the normal. This means that 

5^ - i?o(t) =6 n[i^o(t) ] 
or 


I « 

X - X (t) = 6 

o 


Figure 2 is useful because many of the quantities on it are known. For 
example, it is known that r = c(t - t) since the distance a signal travels is 
sinply the speed multiplied by the time of propagation. In addition, the distance 
the panel moves in this time period is |x^(x) - | which is v(t - x). Thus, 

rM = v(t - x). 

The purpose of this section is to integrate Kg over the area of this arbitrary 
panel. By assuming a planar panel, one can simplify Kg considerably. For example, 
it is obvious that 


cos 0 = M + “ 
n r 


where = v • n/c • Thus, Kg can be written as 


= 


p c^M (M - M^) + p[(1 - M^) (M + 6/r) - (1 - M )M ] 
o n r n r n 


2 3 

r^d - 


ret 


which can be simplified to 




2 2 
M ( p C + p) (M - M ) 
n ^o ^ r 

2 3 




3 3 


ret 
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Therefore 


4ti p(x ,t) 
o 


equation (11) becomes 

= I K as + /j c3S 

f4 ^ fio = 


.Xlm r' C 

6-K) •'-£ J-e 


“2 2 
M (p C + p) (M - M ) 
n *^o r 

2 3 

r (1 - M ) 
r 




r^d - 


dy^ + 0(e) 


ret 


Now assuming that p and v are constant over the e-region yields 
4n p(: 


as + i 


K dS + A K„ dS 
f=0 ^ f=0 ® 


t M (p c + p) lim 
n *^o 

6-^ 


a. r f 

>0 ^-e J- 


e 

e 


M - M 
r 


2 3 

r^1 - M^)'' 


dYi ^2 


ret 


t p lin, f' f" 

6-^0 J-e *^-e 


p^6 


3 3 

r (1 - M ) 
r 


dy^ ^2 + o(e) 


ret 


( 12 ) 


All that is required now to integrate equation (12) is to write the integrands 
in the last two integrals in terras of and y«. From the geometry shown in fig- 

ure 2 and algebra, the propagation vector r can be written in terms of its coitpon- 
ents in the local coordinate system: 


r = (R* sin y\ -R* cos y’/ ^ cos 0) 


where R* , y*, and 0 are defined in figure 2- Using the relations. 


R* sin y* = R sin y 
R* cos y* = R cos y - rM^ 
r cos 0 = 6 t rM 

n 

which are obtained from geometry, one gets 

r = (R sin y, rM^ - R cos y , 6 + rM ) (13) 

t n 
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Therefore the magnitude of $ is governed by 


where 


and 


P^r^ + 2(RM^ cos y - - 5^ = 0 


= 1 - 


2 2 2 
M = M + M 
n t 


Using the quadratic formula gives 


B^r = ( 6 m - RM cos y) + \ ( 6 m - 

^ n t * N n 


RM^ cos y)^ + P^(R^ +6^) 


(14) 


which can be simplified by taking the scalar product of the Mach number vector 


^ = (0, M^, 


and ? (eq. (13)) to get 


rM = M r - RM cos y + 6 m 
r t ' n 


(15) 


Subtracting r from both sides of equation (15) and then multiplying by -1 gives 


r(1-M)=6r+RM cos v - 6 m 
r t ^ n 


Substituting equation (14) gives 


2 2 2 2 

•(1 - M^) = - RM^ cos y) + p (R + 6 ) 


which can be rewritten as 


o o o o o 

r(1 « M ) = Mb R + R M cos Y - 2RM 6M cos y + 6^6 
r t ' t n ’ ^t 


2,2 
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where 




2 

t 


1 


M 


2 

t 


Changing to Cartesian coordinates, since 


r2 = 


2.2 

+ ^2 


and R cos Y ~ y2' 


r(1 - 


M ) = 
r 


\|p"y^ H- p2y2 _ 2M^M^6y^ + p2g2 


( 16 ) 


where 


P 


2 

n 


1 


M 


This is now in a form that allows the integration of equation (12). It remains to 
write in terms of and from equation (15), 


„ ■ - V t 

r r 


(17) 


Substituting equations (16) and (17) in equation (12) gives 

4u p(x ,t) = I K_ dS + /f K dS 
° f=0 ^ f=0 ® 

(6m - My) dy 


n e ^^1 

2 “^ 

e fp y. + 3 - 26 m m. y„ + 


+ M ( p c t p) lim 
o 


2.2'|3/2 


(p"y, + - 26M_^M^y^ + 


i .. (•' r ^ 

■>0 J-E J-e f 6 y . + 


+ p p lim 
6->-0 


6 


( pY '; p ^ y ^ - 26M„H^yj H. p ^ 5 ") 


2.2', 3/2 


+ 0(e) 


which can be simplified to 
4n p( 


= f \ ^ / 


K dS + /f K„ dS 

fio ^ A = 


,P„y^ . p^P, 


5 <iy.| dy^ 




2.2',3/2 


p) lim r f — 2”-2 

6-)-0 J-£ J-e fs y. + 


- M ( p c + 
n t *^0 


^2 *^^1 ^2 


(P^y^ + P^y2 - 26M^M^y2 + P^6^) 


2;,2',3/2 


+ 0(e) 
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Now, to integrate these the following transformations are used: 


I ■ 

5 - ?p„y, 

Because the integrands are symmetric in one gets 


where 


4% p(x ,t) 
o 


* ./ \ ® / 


K as + /f K as 
« R s 

f=0 f=0 


+ (p V + p) lim — ^ 

6^0 p 


~ M M (p c + p) lim 
6^0 


§ 5 an ac 

I *^0 

2 f ari as 

PP^ *^0 Jr\^ (n^ + + p^6^) 


3/2 


2 

T) = -B e - M M^6 
'1 *^n n t 

2 

= 6 e - M M. 6 
'2 n t 

= PP e 

2 '^'^n 


Integrating these in the T)-airection gives (ref. 52, p. 86) 

4n p(x ,t) = lK_as+/JKas + -^(p + p) 

° f=0 ^ f=0 ® o n 

• rk 


X lim 
6->0 


r-T- 

•'o (s + 


612 ^ 


22,2 2 22 

(S + P 6 )\n + S + P 6 


^c, 6 ti as 

I 2 1 

j 222 r"2 2 

•'o (C + P 6 )\ri, + C + 




„2,2 
P 5 


2M M 

n t , 2 

+ (p c + p) 


PP 


2 ‘ o 


n 


X lim 
6^0 


J . _ . (■' 

>0 + P^6^ *^o 




A 2 2 2 2 

+ S + P 6 


+ 0(e) 


0(e) 
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Integrating these in the ^-direction gives (ref. 52, pp- 60 and 89) 


4% p(x ,t) 
o 


I K dS + ^ K dS + — — ( p V + pi 

fio ^ lio = ° “ 

6ti / P e \ 5 ti 

2 -1/ 2 n \ 1 

“ ,| 2 2 2 2 2 2 - TaTT 

6“>’0 2 \ 4*ppe+p6y 1 


X tan 


-i/ ^ 


2 „2„2 2 „2,2 

T) + p p e + p 6 
1 n / 


2M M 

n t, 2 X 

+ T-(p c + p) 

_ 2 O 


/ 12 2 2 2 2 2 \ A |2 2 2\ 

+ p p . + p 6 H- pp - P a ^ 

X lim In , - In . + 0 (e) 

6^0 + p^p^e^ + p^6^ + pp^e/ + p^6y 


where the following transformation is used: 


2 ^2 „2-2 
"’1,2 ^ ^ ^ P ^ 


where ^ signifies t)^ or TI 2 ' whichever is appropriate. Also it is assumed 

that M < 1, so that p > 0. Notice that when c equation (18) reduces to the 

equation of Bisplinghoff et al. (ref. 44, p. 212) for inconpressible flow. 

Now since 

2 

lim n = -lim ri = P e 

6^0 ^ 6^0 " "" 

the logarithmic terms both become zero. Also the inverse tangent terms are equal in 
the limit as 6 0 and behave like 
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which is %/2 in the limit. Thus equation (18) becomes 
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where the integrands are repeated for convenience: 




pcMM +pQ(1-M)+pM cos 9 1 
^o nr r r r ' 

cr( 1 “ M . 

r J ret 


= 


2 2 2 
pcM(M -M)+p[(1-M) cos 9 - (1 - M )M ]l 
^o n r r n 

2 3 

r (1 - M ) 
r 


ret 


Equation (20) represents the governing equation, amenable now to numerical techni- 
ques, for the pressure on a bocty in cortpressible subsonic motion. Tbe theory is 
linearized so that it is e5q>ected to be more and more accurate for thin or slender 
bodies, including most bodies of aeroc^namic interest. For conpleteness, note that 
for the observer inside the hody, sgn(6) = -1; and outside of the body, sgn(6) = 1. 
Thus, equation (2 0) applies to the outside of the body. 

As mentioned in the Introduction (section I. A), these results do not agree with 
those of Morino (ref. 17, appendix C), although a direct cortparison is difficult 
since he uses a velocity-potential formulation. Nevertheless, a few discrepancies 
between the results can be noted. In this section we have calculated the contribu- 
tion from a square e-region. These results are shown to be Mach number dependent 
since the coefficients in equation (20) include Morino claims that the contri- 

bution from the e-region is independent of Mach number, simply 2n, which is the 
incortpressible result. In addition, his e-region is circular. Since the e-regions 
are analogous to the panels in a panel method, it is difficult to understand how one 
can model a body using circular panels. Furthemmore, Farassat (ref. 7) has calcu- 
lated the contribution from a circular e-region for the FW-H equation, and the result 
differs from the result in this section and also from the results of Morino. Thus 
one cannot assume that differently shaped e-regions yield the same result. However, 
the resxilts presented here can be shown to be unchanged for any quadrilateral that is 
symmetric about y 2 / even if one side has zero length. 


D. Keduction to Incoitpressible Aerodynamics 

Since equation (20), especially the integrands and Kg, is relatively com- 

plicated, it is useful to reduce it to its incorrpressible form. If the limit is 
taken as c > all the Mach number terms (M^, M^, M, and l5^) approach zero. 

Furthermore, all the p terms (P, p^, and p^^) approach unity. Thus, = 0 and 
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Therefore, equation (2 0) becomes in the incoitpressible limit 
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This could have been obtained by taking the limit as c of the FW-H equation and 

then solving it using the Green's function of potential theory* 

The first term in the integrand is due to the thickness term of the FW-H equa- 
tion. It is also equivalent to the thickness term of inconpressible aeroc^namics , 
equation (5-85) of reference 44. 

The second term in the integrand is called the "loading** term in aeroacoustics. 
It is directly related to the lifting effects in aerodynamics, for example, equa- 
tion (5-94) of reference 44. 

Throughout this study it has always been useful to refer to the above equation 
for guidance. The reader is encouraged to do so also when something related to the 
full equation is unclear. Of course retarded time effects are not included in the 
above equation because of the infinite speed of sound propagation. 


III. COMPUTATIONAL METHOD 
A. Approximating Body by a Finite Number of Panels 

The cortplexity of equation (2 0) precludes an analytic solution for most, if not 
all, bodies of interest. Therefore one must use approximate numerical methods. The 
analysis in the previous section showed that the solution can be approximated to 
0( e) where is the area of a hole removed from the region of integration. This 

procedure was shown to be analogous to taking the Cauchy principal value which is 
also an approximation that is accurate to the order of the size of the hole deleted 
from the integration. 

In the past, to predict the noise due to bodies in motion, one needed the sur- 
face pressure. Now, using equation (20), one can calculate the pressure on a small 
region of the body, given the pressure on the remainder of the body. By discretizing 
the surface, one can develop a system of algebraic equations whose solution gives the 
pressure everywhere on the body. To do this, the body must be approximated by a 
finite nximber of planar elements, or panels. The pressure is assumed constant over 
each element. The solution approaches the exact solution as e 0. 

As mentioned earlier, the e-region corresponds to a small area of the surface 
immediately around the observer. To develop the system of algebraic equations, the 
observer must be located on each panel of the body successively. Each location of 
the observer yields an algebraic equation. 

Recognizing the above, one can proceed in several ways to solve the equation 
numerically. They differ in how the dependent variable is assumed to behave. One 
method that has been used on similar equations is a collocation procedure (ref. 53) 
using global shape or loading functions. This method is common in lifting surface 
theory (refs. 54 and 55) , where the loading functions are chosen to exactly satisfy 
the leading-edge, trai ling-edge, and tip conditions and to produce satisfactory 
numerical results. Once these functions are selected, there are still integrations 
to perform, although these are now over known functions and may permit analytic 
evaluation. Otherwise quadrature is necessary. 

The collocation method was not selected for the present work because of its lack 
of generality; the shape functions required depend on the boc^ shape. For example. 
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one may have to use a different set of functions for a fuselage, a propeller, and a 
wing - not to mention for more complicated shapes such as highly swept propeller 
blades. 

Since the purpose of the numerical procedure in the present work was to verify 
the feasibility of the method and not to develop a production computer program, it 
seemed appropriate to solve the equation directly rather than to introduce additional 
complications. Therefore, the computer program that was developed uses only quadra- 
ture. The dependent variable is not represented by global shape functions. Of 
course, more efficient or sophisticated methods may be available for particular 
applications. As an initial solution method however, this approach brought out the 
siibtleties and pitfalls of the equation better than any other. A disadvantage of the 
method is that it produces a large system of equations, but they are strongly diag- 
onal. Therefore the system of equations is readily solved. 

Quadrature is used over each element. As mentioned earlier, the observer, at 
x^, is positioned on each element successively (at time t). The observer is 
stationary, so one must specify the time on the body. For the observer on a given 
panel, the numerical integrations are carried out over all the other panels. The 
theoretical developments presented earlier have alrea<^ calculated the effect that 
the observer panel has on itself. 

For a body approximated by N panels as in figure 3, the approximate form of 
equation (20) is 



Figure 3.- Bo<^ approximated by panels. 


where 
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and subscripts i and j signify the observer and source panels, respectively; 
the observer and source are assumed to be at the centroid of the respective panels. 
Equation (21) is one equation for N unknovms, and p^ (j = 1/ •••, N, 

j ^ i) . Letting i take on the values from 1 to N gives a system of N equa- 
tions for N unknowns. 

This system of equations can be e5q>ressed by a matrix multiplying a vector of 
unknowns equal to another vector: The coefficients of p^ are the diagonal terms of 

the matrix. The integrals over Kp are the off-diagonal terms. The right hand side 
of equation (21) yields a unique value for index i. 


B. ^adrature Formulas 

Legendr e-Gauss quadrature formulas are used to calculate the integrals 
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K dS 
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Portions of each integral can be ill-behaved because of the 1/r terms in it. Even 
though the effect an element has on itself has been calculated, two different ele- 
ments may be very close together, for example, on the upper and lower surfaces of a 
thin airfoil, especially at the trailing edge. One -dimensional Legendre -Gauss quad- 
rature is exact for polynomials of order (2n - 1) where n is the number of nodes. 
For the above integrals this may not be sufficient. They may require very high order 
polynomials to approximate their behavior because of the negative exponent on r. 

Just as for the analytical part of this work, the numerical solution would be 
trivial if not for the ill-behaved nature of the integrands. Since this does repre- 
sent the most difficult aspect of the computational method, two simple examples are 
given to illustrate the singular behavior of the integrals. 

The ill-behaved terms in both and become simply 

pressible flow. Thus the integrals become proportional to 


cos 9/r^ for incom- 
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cos 0 dS 

2 
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Now for these examiples, assume that the integration region (source panel) is a unit 
square and the observer is located along the normal directly above the center of the 
panel, as in figure 4. This is approximately the orientation of the two panels that 
lie directly across from one another on the upper and lower surfaces of an airfoil. 
The observer is at x and the source point is at y. Iherefore, I represents the 
effect that the unit square has on the observer. 
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pressible flow. 


Note that 


cos 9 

2 

r 


(?? - 4 ^ 


and the integral becomes 


1 ( 6 ) 




_^2 

t 62)3/2 


( 22 ) 


which can be evaluated exactly as 


1 ( 6 ) 


4 tan 



(23) 


This integral illustrates the difficulty in solving equation (21) using quadrature. 
Notice that the integration is through the point which the integrand 

behaves like 1/6^. For small 6, it is very difficult to approximate with a poly- 
nomial a function that grows this rapidly. Were it not for this behavior for two 
panels close together, equation (21) would be relatively easy to solve. Also, 
because bodies of aerot^namic interest are usually thin, panels do lie close 
together. 


Two approximations of the integral 1(6) using Legendre-Gauss quadrature 
(ref. 56) and different numbers of nodes are now presented. The one-dimensional 
formulas over a surface give 


1 ( 6 ) 


n* n* 
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where n’ is the number of nodes in each direction, and Hj are weight coeffi- 
cients, and (C^ i'^2 node locations (given). The node locations and 

weight coefficients are given in appendix A- When n* =2, the above formula simply 
becomes 


1 ( 6 ) 


46 

[2(0.57735)^ + 


(24) 


where all the weight coefficients are unity and all the nodes are at ±0.577350. 


The inportant question is how accurately equation (24) models the exact solu- 
tion, equation (23), for the various values of 6. One would expect it to be inaccu- 
rate for very small 6 because the integrand varies so rapidly. 


Figure 5 is a graph conparing the exact value of the integral 1(6) and numeri- 
cal approximations to it. The abscissa uses the relation 1/6^ because this is pro- 
portional to 
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ret 
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Figure 5.- Coitparison of Legendre -Gauss quadrature 
and exact integration for sample integral. 
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which is used in the conputer program to determine how many quadrature points to 
use. In addition to equations (23) and ( 24 ), the results from the use of an 
8 -point (64 points over a surface) quadrature formula are presented. For values of 
1/6^ < 0.19, the 2-point formulas are accurate to at least 1 percent. This same 
degree of accuracy can be obtained for even larger values of 1 / 6 ^ using the 8 -point 
formula. These approximate results diverge very rapidly from the exact value for 5 
near zero. 

These corparisons not only illustrate the behavior of the integrals but also 
provide a means to determine what order quadrature should be used (i.e., how many 
nodes) , given the area of a panel and the distance between the panel and the 
observer. To achieve reasonable efficiency and accuracy, several orders of quadra- 
ture formulas are used in the current conputer program, which actually determines how 
many nodes to use for each integration. The maximum is 64 nodes over a panel. The 
minimum, for panels and observers far apart, is a simple 1 -point rectangular rule. 
Further studies were conducted using quadrature formulas that were exact for trigono- 
metric functions, but they offered no significant improvement in accuracy or 
efficiency. 

In the conputer program, for panels that have Fatio > 16, an exact form of the 
solution is used (eq. (18)), that is, the form of the equation before the limit as 
6 > 0 was taken. But because the observer is assumed (in the development of 
eq. (18)) to lie along the normal to the source panel, this formula can be used only 
for panels oriented in this manner. This does not represent a problem though, since 
the panels on an aerodynamic body that are closest together (e.g. , at the trailing 
edge) are so oriented. 


C. Jacobians and Mapping of Elements 

Since Legendre -Gauss quadrature requires the integration region to extend from 
-1 to 1, each element must be mapped to a unit square, by using standard finite- 
element transformations. Of course, to integrate over the mapped element, one must 
also calculate the Jacobian of the transformation. This section briefly describes 
the methods used to perform these operations. 

The elements that make up the bodies are all in three-dimensional space. They 
are in motion, but that does not enter into the following considerations. The 
retarded time aspects affect only the integrands and not the integration surface 
directly. 

The required integrals are of the form 


^ = J i K(x,y,z) dS 
Panel 


To reduce confusion, the Cartesian coordinate system is denoted by x, y, and z 
rather than by x^ , X 2 , and x^. If this panel is mapped to the ^ 13 ”^ plane of 

the rj-frame, it becomes (ref. 57) 
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where 
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The panels are mapped to the ^-frame with linear finite element shape functions 
(ref. 58) of the form. 
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where and (i = 1 to 4) represent the vector coirponents of the corners 

of the source panels in the original coordinate system. Further, 
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The Jacobian at any particular value of 
ating the above formulas and using equation (2 5) . 


and r \2 is obtained by differenti- 
For example. 
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These formulas are used in the subroutines JACOBI and NODE. A flowchart of all the 
subroutines in the computer program is given in appendix B. In JACOBI, the Jacobian 
at every node on every element is calculated and stored. Since it is not known a 
priori how many nodes are required to perform the integrations, the Jacobians for 
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each element are calculated for 1-, 4-, 16-, and 64-node quadrature* Thus, the 
Jacobian is calculated for 85 points on each panel. All these Jacobians are stored 
except the 64-node values, which are calculated when needed. This means that 21N 
Jacobians are stored, where N is the number of elements. 

Given a node and subroutine NODE calculates x, y, and z from equa- 

tions (26). Values of r]^ and r ]2 corresponding to Legendr e-Gauss quadrature (see 
appendix A) are stored in DATA statements at the beginning of the main program. 

It is important to point out that the Jacobians must be calculated only once for 
a given body. Changes in angle of attack or velocity do not affect them. 


D. Retarded Time Calculation 

Retarded time is another quantity that must be calculated numerically. It is 
typically governed by a transcendental equation that is difficult to solve analyti- 
cally. In this section, the governing equation is derived for retarded time for a 
body that is moving rectilinear ly (along the z-axis) and spinning (about the z-axis) 
(see fig. 6) . This motion is not the most general type possible but is adequate to 
test the theory for several types of bodies. For instance, because the spinning is 
about the same axis as the rectilinear motion, the motion of a helicopter blade in 
forward flight could not be represented with the governing equation derived here. 
However, other types of motions should be relatively easy to include in future 
programs. 

The quantity to be calculated is the retarded time distance between an 

observer and a source point. This is the r that appears throughout this study, 
for example, in equation (8). It represents the propagation distance for a signal 
emitted from a source that is in motion. The observer is stationary, as has been 
assumed throughout the report. 
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Figure 6.- Velocities and position vectors used 
in retarded time calculations. 



In Cartesian coordinates, the observer is located at 


o o o o 


and the source in motion is located at 


r(T) = [x cosCcot) - y sin(wT), x sin(a}T) + y cos(a)T), z + Ut] 
s s s s s s 

where t = t - r/c, o) is the angular velocity of the body, and U is the recti- 
linear velocity of the bo<^. When t = 0, the source is at [x , y , z ]. 

S o S 

Now from its definition. 
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Since only stea<^ motion is being considered, the choice of a particular value of t 
is arbitrary. The pressure is not changing with time (in these body-fixed coordi- 
nates) • Therefore, setting t = 0 gives 
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This is the transcendental equation governing subroutine RADIUS (see 

appendix B) this equation is solved with a Newton-Raphson technique (ref. 56) . 
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Once determined, it can be used to calculate the vector quantity 

from 



E. Coiqjuter Program 

This section describes the conplete cortputer program. The major tasks of the 
program are described in the order in which they are performed, as shown in the flow- 
chart in figures 7 and 8. 



Figure 7.- Conputer program 
flowchart, part I. 
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The first major task of the program is to read an input file that describes the 
body and its velocity ( box 1 ) . The body is described by a finite number of coordi- 
nate points in three-dimensional space. These must be input properly to form the 
elements. They are input in the order shovm in figure 9 , which shows a wedge-shaped 
airfoil. The corresponding element numbers are shown in figure 10. Each cross 
section must have the same number of points, although this would not have to be so 
in future programs. Note that to "close" each cross section, the first and last 


z 



Figure 9.- Proper order of input 
coordinate points. 


z 

i 
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Figure 10.- Numbering of elements. 


36 



coordinate points are coincident. A good way to check the input points is to examine 
the normal vectors that the program calculates. Quite often if the points are input 
in an improper manner, the normals point into the body. 

The next step in the code is to transform the coordinates to account for the 
angle of attack ( box 2 ) • This could also be used to change a set of coordinates 
describing a rotating blade to those describing a wing in rectilinear motion 
(a = 90® ) • Since rectilinear motion is always in the direction of the z-axis and 
spinning occurs about this axis, the hody would have to be rotated to switch roles 
from a blade to a wing. 

The current program is designed to run only one angle of attack at a time. 
Ideally, for a production code, another loop would be added to calculate multi- 
ple right hand sides (RHS) of equation (21). Remember it is necessary to calcu- 
late the left hand side only once for a given Mach number, but the RHS varies for 
each angle of attack. This procedure would be especially efficient since most 
linear equation solvers are equipped to handle multiple RHS and give multiple solu- 
tions. Ihe matrix would be manipulated only once. 

The next major section of the program ( box 3 , fig. 7) is a loop that calculates 
normal vectors, Jacobians, and centroids of each element. The normal vectors are 
calculated in the main program using a technique described by Hess (ref. 27), which 
simply takes the vector cross product of the diagonals of each panel. 

Jacobians are calculated in the subroutine called JACOBI (see appendix B) using 
the formulas presented in section III.C. The centroids are calculated in subroutine 
CENTER by integrating over each panel. Since the nodes and Jacobians are already 
calculated and stored, this is relatively straightforward. The area of any panel 
can be calculated simply by adding the four Jacobians calculated for the two-point 
be gen dr e-Gauss quadrature for that element, since the weight functions for the two- 
point formulas are unity. 

At this point in the program all the preliminaries have been conpleted. The 
remainder of the program is shown in figure 8, the second part of the flowchart. The 
first DO-loop, over i ( box 4 ), places the observer on each element successively. 

For each location of the observer, all the other panels are successively treated as 
sources by the next DO-loop over j ( box 5 ) . Quadrature is performed over each 
source panel. 

Before the quadrature is performed, the program must determine how many quadra- 
ture points to use. For each value of i and j, the quantity 
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is calculated, where ^ and ^ are taken as the centroids of the i and j 
panels, respectively. As shown earlier in section III.B, the mamber of nodes 
required for accurate quadrature is highly dependent on this quantity. Iherefore, 
the program calculates this quantity in subroutine TEST and determines how many 
quadrature points to use. The area of the source panel is used for dS, of course. 


Once the appropriate number of nodes to use is determined, the actual quadrature 
calculation begins at the DO-loop over k (box 6). The number of nodes it will use 
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for that particular integration is Subroutine NODE determines from equa- 

tion (26) the location of the node in boc^ coordinates from the Gaussian nodal coord- 
inates. Then subroutine RADIUS calculates the retarded time propagation distance 
^ref Subroutine MACH calculates M, M^, and M^. With these quantities cal- 
culated, the integrands and can be conputed for equation (21). By multi- 

plying the integrands, the weight coefficients, and the Jacobian at each node, and 
then adding these products from each node, the integration is complete over a parti- 
cular panel. 

After the DO- loop over i is cortpleted, the system of algebraic equations is 
completely developed. The solution to the system of equations is obtained using 
subroutine INVERT, which, in addition, conditions the equations to satisfy a Kutta 
condition. Ihis conditioning is described in the next section. The final step of 
solving these equations is performed using a NASA Langley Research Center library 
subroutine called GELIM (ref. 59), which uses an upper or lower decomposition as 
described in reference 56. 


IV. LIFTING BODIES 

Determining the surface pressure on lifting bodies is, of course, much more 
difficult than on nonlifting bodies because lift is an inherently viscous effect. To 
approximate viscous effects with an inviscid theory, one must obviously use not only 
physics but also e>q>erimental results to extend the theory. But, as Hess (ref. 26) 
states: 


It is important to realize that any such formulation is simply 
a. . . .model of real lifting flow, and that the two flows are not 
necessarily related in any fundamental way. 

In this section the major differences between inviscid and viscous theory are briefly 
described, as they relate to aeroc^namics. After that a method for using the present 
theory to approximate the effects due to viscosity is presented. These discussions 
are more qualitative than quantitative. 

The linearized, inviscid equations used herein neglect several important fea- 
tures of the actual flow. First, the flow is assumed to slip over the body, which 
means that no tangential velocity boundary condition is prescribed. This immediately 
leads to a question of uniqueness, for there are an infinite number of velocity 
fields that satisfy the normal velocity boundary condition. The Kutta, or trailing 
edge, condition is usually used to determine which of these velocity fields is appro- 
priate, since it is known that viscosity prevents the flow from turning the sharp 
angle at the trailing edge. 

Neglecting the no-slip condition also means that there would be no boundary 
layer or wake in this inviscid flow, although there could be a layer of zero thick- 
ness over which a tangential velocity jvmp occurs. Since the pressure is continuous 
across such a layer, this jump does not affect the formulation when using the pres- 
sure as the dependent variable, as mentioned in section I.A. 

A convenient, and corrputationally inexpensive, way to account for the boundary 
layer occurring in real flows is to use the displacement effect (refs. 35, 60, 
and 61) via an iterative technique. First, the inviscid calculations are performed 
with no boundary layer assumed. Then, from these results, the displacement thick- 
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ness is calculated. rflie boundary layer and its effect on the streamlines is then 
accounted for by enlarging the actual body by an amount equal to the displacement 
thickness. 

Another difficulty in using the FW-H equation to predict lift is that there is 
no built-in mechanism to generate lift. The circulation is assumed to be zero. The 
well-known panel methods distribute line vortices on or inside the body/ where their 
strengths are determined by satisfying a Kutta condition. These vortices simply 
provide a means of altering the onset flow by providing a circulatory contribution to 
the flow. Of course, one must still satisfy the normal velocity boundary condition. 
Using the FW-H equation as the governing equation means that there is no simple way 
to include line vortices. Fortunately, there is a way to condition the final system 
of linear equations and obtain the effects due to lift without using the concept of 
vortices distributed within the body. 

To understand fully the technique used to condition the equations, one must con- 
sider what the flow field is really like around a lifting airfoil. A major differ- 
ence between a lifting and a nonlifting airfoil is the location of the leading^edge 
stagnation point. Conversely, the Kutta condition tells us that the trailing-edge 
stagnation point is in the same location for lifting and nonlifting airfoils. Fur- 
thermore, even the streamlines near the trailing edge are remarkably similar for 
lifting and nonlifting airfoils. Thus, the trailing edge experiences essentially the 
same flow whether the body is at an angle of attack or not. Van Holten (ref. 12) 
essentially implies this fact when he advocates placing only the leading edge at an 
angle of attack. He uses the acceleration potential and finds that using this 
"variable- geometry** concept forces his inviscid results into close agreement with 
real lifting flows. Apparently, these ideas can be traced back to Lanchester. 

To illustrate why Van Holten* s suggestion works, several numerical examples are 
presented for a simple body made up of only six elements. By presenting numerical 
results for a simple body f the procedures should become clearer. Thus, equation (20) 
is solved numerically by the computer program described in the previous section for 
the body shown in figure 11, a crude six-element model of a wing with an aspect ratio 
of 3. The profile, which is roughly that of an NACA 0012 airfoil, is shown in fig- 
ure 12. Using this simple body allows one to investigate the numerical procedure and 
to illustrate how the effects due to lift are included. The aspect ratio is deliber- 
ately kept small so that only a few elements need be used. 



Figure 11.- Simple six-element airfoil with 
element numbers. 
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Figure 12.- Profile of simple six-element 
airfoil. 


'The first example is for the body in very low-speed rectilinear motion 
(M = 0.03) at an angle of attack of zero. The resulting system of equations is 


2n 

0.19 

0.07 

0.07 

0.43 

5.43 


c 


— 0.66 

0.20 

2n 

.19 

.25 

4.92 

.61 


c 

P2 


-1.70 

.17 

.28 

2n 

.10 

.82 

.22 


c 

^3 


6.17 

.22 

.82 

.10 

2n 

^28 

%17 


c 

P 4 


6.17 

.61 

4.92 

.25 

.19 

2n 

.20 


c 

P5 


-1.70 

5.43 

.43 

.07 

.07 

.19 

2n 


c 

_^6 


— .66 


where 


is the pressure coefficient on ith element. The solution to equa- 


tions (27) is shown in figure 13. It differs markedly from the results for the two- 
dimensional airfoil (from ref. 62) which are also shown, but that is to be ejected 
since the aspect ratio is only 3 and the paneling is very crude. However, it does 
show the proper trend, which is all that is desired since these discussions are 
designed only to illustrate the general behavior of the equations. An alternative 
approach of discussing these ideas using the equations themselves is possible but 
would probably not be as successful in instilling an intuitive feeling since the 
discussion would be very abstract. Much can be learned concerning equation (20) 
from the examples included herein. 


Getting back to the exairple, notice that in equations (27) , there is a symmetry 
defined by: 


a. . 
ID 


^7-i,7-j 


b. 

1 


b 


7-i 


where i and j range from 1 to 6. This means that the pressure on the upper and 
lower surfaces is governed by the same equation, since the airfoil is symmetric and 
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at zero angle of attack, so that the upper and lower surfaces experience the same 
flow. It also means that 


c = c 


Pressure 

coefficient 

c 


P 



0 .2 .4 .6 .8 1.0 


Distance along chord, x/c 


Figure 13.- Pressure distribution on simple six-element 

airfoil. 


Another irrportant point is that some of the off-diagonal terms in the matrix are 
nearly as large as the diagonal terras, especially for the panels at the trailing edge 
(the first and last rows of the matrix) , because they are so close together. The 
effect one of these panels has on another is nearly the same as the effect one has 
on itself. In addition, these large off-diagonal terms are the most difficult to 
calculate, as discussed in section III.B. For an infinitely thin bo<^, these large 
off-diagonal terms would be equal to the diagonal terms. Kuo and Morino (ref. 63) 
observed this fact also, but found that it did not cause any serious numerical prob- 
lems for bodies of interest. It is also of interest to note that for an infinitely 
thin bo<^ the right-hand- side vector would be zero. Therefore, one would obtain the 
result 


c = -c 

^upper ^lower 


and consequently zero lift. 
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The really interesting features appear when the hody is placed at an angle of 
attack. Using the same airfoil as in the previous example but at an angle of 10® 
yields the following system of equations: 


271 

0.19 

0.07 

0.07 

0.43 

5.43 


c 


0.24 

0.20 

2ti: 

.19 

.25 

4.92 

.61 


c 

^2 


-2.01 

.17 

.28 

2% 

.10 

.82 

.22 


c 

^3 


3.85 

.22 

.82 

.10 

2tc 

.28 

.17 


c 

P 4 


8.31 

.61 

4.92 

.25 

.19 

2'ji 

.20 


c 

P 5 


-1.17 

5.43 

.43 

.07 

.07 

.19 

2ti 


c 

Pe 


-1.47 















Notice that the matrix is unchanged from the a = 0® matrix, since its elements 
depend only on the relation between the various points on the bo^. Of course, for 
different Mach numbers the matrix changes, but it is not angle -of -attack dependent. 


The angle of attack does significantly alter the RHS vector of the system of 
equations. This change results in the solution presented in figure 14. The large 
jump in pressure across the trailing edge is representative of flows over airfoils 


Pressure 

coefficient 
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P 



Figure 14.- Pressure distribution on simple six-element 
airfoil when a = 10°. 


before applying a Kutta condition. Satisfying the Kutta condition means forcing the 
pressure on the upper and lower surfaces to be equal at the trailing edge. To under- 
stand how to accomplish this, one must remember that the Kutta condition, as stated 
earlier, forces the flow at the trailing edge to be as it was when the airfoil was 
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at a = 0® , at which the conditions on the upper and lower surfaces were the same. 
The only difference in the system of equations though is the RHS vector. When 
a = 0®/ = b , ^2 ” ^5^ ^3 ” ^4^ these no longer hold true. How- 

ever, they do iflustrate a method of applying a Kutta condition. If b^ = b^ 
were enforced for the system, then the trailing edge of the airfoil would ejqoerience 
the same flow on the upper and lower surfaces. Tiierefore, in equations (28) replace 
b^ and b^ by 



which gives 
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— — 


The solution of this system is shown in figure 15. By using averaged values for b^^ 
near the trailing edge, a Kutta-type condition was enforced which produces a net 
force on the body due to the difference in pressure on the upper and lower surfaces. 


Pressure 

coefficient 
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Figure 15.- Pressure distribution on simple six-element 
airfoil when a == 10® with Kutta condition. 
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One aspect of this conditioning technique that was not illustrated very well by 
this simple example is that it must be applied over a large portion of the chord. 

In this exartple the trailing-edge panel was very large, covering 50 percent of the 
chord. If many elements had been used along the chord, this conditioning would have 
had to be applied to several, unlike velocity-potential panel methods where vorticity 
is distributed inside the body. When using the vorticity method, one must apply a 
Kutta condition only on the panels adjacent to the trailing edge because the vorti- 
city tends to couple strongly the leading- and trailing-edge algebraic equations 
(in fact the equations for the whole airfoil). When using the pressure formulation, 
there is no strong coupling between the leading and trailing edges. In fact, the 
panels on the upper and lower surfaces of the trailing edge affect each other so 
strongly that they are almost unaffected even by adjacent panels. Thus, applying the 
Kutta condition only to the single pair of panels at the trailing edge would have 
almost no effect on the next set of panels forward from the trailing edge. There- 
fore, the conditioning must be applied to a large section on the aft of the airfoil. 
From preliminary studies on uncambered, symmetric airfoils, it has been determined 
that this conditioning of the RHS vector must be applied to all panels from the 
trailing edge to the thickest part of the airfoil. Thus, for an NACA 0012 airfoil, 
with many panels along the chord, one would condition the vector over 70 percent of 
the chord, since this airfoil is thickest at 70 percent of the chord (from the trail- 
ing edge). As an exaitple, if one modeled the airfoil by 10 panels on the upper sur- 
face and 10 on the lower, where each panel covered 10 percent of the chord, the 
vector would be conditioned as follows: 


= ^20 = 


1 20 


^2 = ^19 = 


2 19 


K = Ka = 

7 14 


t> + b 
7 14 


The above technique has an effect coitpletely analogous to that of using Van Holten's 
suggestion, except that its implementation is much simpler because the body is not 
considered to be deformed in any manner. 

It is also interesting to note that using the displacement thickness effect also 
conditions the equations in a similar manner. By thickening the body by an amount 
equal to the displacement thickness, one moves the two trailing-edge panels (upper 
and lower surfaces) farther apart. This tends to move the values of b^ on the RHS 
of equations (28) closer together. The larger the displacement thickness, the closer 
the values become. 
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V. RESULTS 


This section presents mjitierical results for several different types of bodies in 
various combinations of rectilinear and angular motion. Using the cortputer program 
described in section III, the surface pressure was calculated for prolate ellipsoids, 
wings, and rotating blades. These different types of bodies were used to illustrate 
the generality of the method. 


A. Prolate Ellipsoids 

Prolate ellipsoids with fineness ratios of 0.05, 0.10, and 0.25 were modeled 
in low-speed motion. These were chosen because of the availability of theoretical 
(inconpressible) and ej<perimental results with which to conpare (ref. 64). 

First, the results for a 5-percent-thick ellipsoid are presented. The paneling 
used is illustrated in figure 16. Hiere were 264 panels and symmetry was not 
ejqploited. If one takes advantage of the numerous planes of symmetry, the number of 
computations and the amount of storage can be reduced substantially. Panels were 
concentrated near the ends of the ellipsoid because of the extremely rapid change in 
pressure occurring there. Remember that the pressure is assumed constant over each 
panel. 


The numerical and experimental results are cortpared in figure 17, which is a 
graph of pressure coefficient versus axial distance. The solid line represents 
e:jq)erimental results (ref. 64), and the symbols represent the current numerical 
results. Agreement is generally good except near the stagnation point. Poor agree- 
ment there is expected because of the use of linearized theory. Near the ends of the 
ellipsoid the flow is essentially that for a blunt body, and linearized theory is 
usually valid to within a few percent of the chord from the leading edge (ref. 44, 
p. 209). For more detailed descriptions of the flow near the stagnation point one 
could use the methods of Van Dyke (ref. 65) or include the quadrupole teirms , which 
include the nonlinear effects, in the present formulation. Neither would be trivial 
to incorporate into the current scheme. 

Figure 18 shows the paneling used for a 10-percent-thick ellipsoid, and results 
similar to those just discussed are shown in figure 19. Figures 20 and 21 show the 
same information for a 25-percent-thick ellipsoid. From these figures it is obvious 
that the nonlinear terms become more and more important for the thicker bodies. Of 
course, these are extreme examples in that most aerodynamic shapes have thicknesses 
of 15 percent or less, especially those designed for higher speeds. For a detailed 
discussion of the implications of nonlinearity in the pressure formulation, see sec- 
tion I. Further, ellipsoids are a severe test of linear theory because of their 
blunt ends. Much more accurate results would be expected for shapes such as ogives 
and wings. 
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Figure 16* - Paneling used to model 5-percent-thick 

ellipsoid. 



Distance along axis, z/L 

Figure 17.- Pressure distribution on 5 -per cent- thick 

ellipsoid. 




Figure 20. - Paneling used to model 25-percent-thick 

ellipsoid. 



-,25 ^ 
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Figure 21.- Pressure distribution on 25-percent-thick 

ellipsoid. 


B, Finite Wing 


The next example is a finite wing with an NACA 0008 cross section and an aspect 
ratio of 5. The paneling used and its cross section are shown in figure 22. There 
were 242 panels with no special precautions taken near the tip. In fact, the tips 
are not paneled at all; they are actually open. Including these panels had an insig- 
nificant effect on the panels away from the tips. And since this section conpares 
numerical results with two-dimensional section airfoil theory, only results from the 
center of this finite wing are presented. 

The numerical results from the program described in section III are presented in 
figure 23. Also shown are two-dimensional potential theory results (incompressible) 
(ref. 62). The numerical results were obtained for a Mach number of 0.01 and coirpare 
very well with the theoretical results. The finite aspect ratio may account for the 
numerical results being slightly lower than the theoretical curve. The numerical 
results for M = 0.06 are also included. The pressure coefficients are higher than 
those at M = 0.01. 



Figure 22.- Paneling used to model NACA 0008 wing. 
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2-D incompressible 
flow theory (ref. 62) 

Present method (M = 0.01) 



Figure 23.- Pressure distribution on NACA 0008 wing 



C. Helicopter Rotor Blade 


In this section mutierical results are presented for a rotating NACA 0012 blade. 
This airfoil section is common in helicopter applications. The results are coitpared 
with the experiments of Gray et al. (ref. 66). 

The blade was approximated by 264 finite panels as shown in figures 24 and 25. 
Panel sizes were decreased smoothly from root to tip. Blade length was 6.5 chord 
lengths. Tip Mach number was 0.25 and angle of attack was zero. 



Figure 24.- Planform view of paneling used to model NACA 0012 

helicopter rotor. 



Figures 26, 27, and 28 conpare the variation in the pressure coefficient along 
the chord for three spanwise locations: 94, 98, and 99.5 percent. The solid line 

represents two-dimensional airfoil theory (ref. 62) . The circular symbols are exper- 
imental results (ref. 66) and the squares are the present numerical calculations. 

The leading edge is obviously at x/c = 0. 
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These results agree well with experiment. Notice how close to the two- 
dimensional theory the results are at the 94-percent span location. Three- 
dimensional effects become obvious as one moves closer to the tip. Of course at 
the leading edge the well-known square-root singularity appears. The solution in 
this region can be approximated with the methods described in reference 65. 

Figures 29, 30, and 31 present similar results for the blade at an angle of 
attack. The concepts discussed in section IV allow the calculation of surface pres- 
sures on lifting bodies. As described there, the inhomogeneous terms of the system 
of algebraic equations must be conditioned by replacing the Inhomogeneous terms (for 
the panels on the upper and lower surfaces at the trailing edge) by the average of 
the two values. Tliis must be done for all panels from the trailing edge to the 
thickest part of the cross section. Here too the nvimerical results agree well with 
ejq>eriment (solid lines) . The circular and square symbols represent the numerical 
results from the upper and lower surfaces, respectively. 

Now, surface pressure calculations from an NACA 0006 rotor at various Mach num- 
bers are presented. Figure 32 shows the pressure coefficient at 94 percent span 
versus distance along the chord for tip Mach numbers of 0.5, 0.7, and 0.9. The solid 
curve indicates the two-dimensional section result for incompressible flow. Notice 
that the numerical results are markedly different from those that would be expected 
from two-dimensional compressible flow. The Prandtl-Glauert rule indicates that the 
line should shift upward for increasing Mach number according to 


c 

P 


c 


P 


o 



where Cp^ is the incorrpressible pressure coefficient. But because of the highly 

three-dimensional nature of flow around a rotating blade, this rule does not apply. 
This is verified by the corrputational results, which deviate significantly from what 
the Prandtl-Glauert rule would predict. Note that for large aspect ratio wings in 
rectilinear motion, the cortputational method presented herein does predict behavior 
similar to that of the Prandtl-Glauert rule. 

One final point concerns the amount of cortputation time that was used on 
these calculations. It is relatively small considering the details obtained and 
the complexities involved. The most time-consuming case, the NACA 0006 rotor at 
M^j^p = 0.9, required nearly 40 minutes of CPU time on a CDC CYBER 175 conputer. The 
same blade at ” 0*01 required 21 CPU minutes. Remember that this was merely a 

pilot conputer program, and the efficiency could be enhanced by optimizing the code. 
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Figure 26.- Pressure distribution on NACA 0012 rotor at 94 percent 
span when a = 0*^* M. = 0.2 5* 
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Figure 27.- Pressure distribution on NACA 0012 rotor at 98 percent 
span when a = 0® . M. . = 0.2 5. 
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Figure 28.- Pressure distribution on NACA 0012 rotor at 99.5 percent 

span when a = 0® . M . . - 0.25. 
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Figure 30.- Pressure distribution on NACA 0012 rotor at 98 percent 
span when a = 6.18®. M. = 0.25. 
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Figure 31.- Pressure distribution on NACA 0012 rotor at 99.5 percent 
span when a = 6.18®. ~ 0.25. 





VI. CONCLUSION 


This report describes a method for calculating the surface pressure on arbitrary 
bodies in conpressible flow. It is based on the acoustic analogy approach of Ffowcs 
Williams and Hawkings and the equations due to Farassat. Its major contributions 
have been in interpreting Farassat* s equation when the observation point is located 
on the surface of the body and in solving the resulting equation numerically. 

Since Farassat*s equation ,is based on linearized, inviscid theory, the method 
presented here of course has its limitations. These have been pointed out whenever 
possible in the text. For bodies moving in rectilinear motion there may be more 
appropriate theories in which some nonlinear effects may be included relatively 
easily. For bodies such as propellers and helicopter rotors, for which the problem 
is very complex, the present results are very encouraging. The method is based on 
first principles and is relatively easy to coitprehend, yet it provides an ine>q>ensive 
way to obtain very useful information. In addition, because the Ffowcs Williams- 
Hawkings equation does contain viscous and nonlinear terms, there is hope for includ- 
ing these effects in the future. 

The method described herein contains several ideas that are somewhat novel to 
aerodynamic theory. First, compressibility is accounted for exactly as it occurs in 
nature, in terms of retarded time. Compressibility manifests itself as a finite 
speed of propagation of disturbances, and in the current theory the distance the 
signal actually travels is calculated. In the past this was an insurmountable com- 
plication; but because of recent advancements in computer technology, it is now 
possible to model these effects. Furthermore, the effects due to lift are included 
without recourse to distributions of vorticity. A method is presented for condition- 
ing the final system of algebraic equations (generated by the conputer program) to 
account for lift via a Kutta, or trailing-edge, condition. This conditioning amounts 
to setting the pressure equal on the upper and lower surface of the trailing edge. 
Vorticity is not required because the theory is formulated in terms of the pressure. 


Langley Research Center 

National Aeronautics and Space Administration 
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APPENDIX A 


LEGENDRE-GAUSS QUADRATURE FORMULAS 


Notation: 


Node locations, 

N = 

0.5773502691 8962576451 


0.7745966692 4148337704 
0.0000000000 0000000000 


0.8611363115 9405257523 
0.3399810435 8485626480 


0.9061798459 3866399280 
0.5384693101 0568309104 
0.0000000000 0000000000 


0.9324695142 0315202781 
0.6612093864 6626451366 
0.2386191860 8319690863 

N = 

0.9491079123 4275852453 
0.7415311855 9939443986 
0.4058451513 7739716691 
0.0000000000 0000000000 


0.9602898564 9753623168 
0.7966664774 1362673959 
0.5255324099 1632898582 
0.1834346424 9564980494 


0.9739065285 1717172008 
0.8650633666 8898451073 
0.6794095682 9902440623 
0.4333953941 2924719080 
0.1488743389 8163121088 


0.9815606342 4671925069 
0.9041172563 7047485668 
0.7699026741 9430468704 
0.5873179542 8661744730 
0.3678314989 9818019375 
0.1252334085 1146891547 


Weight coefficients, 

2 

(1)0.1 000000000 0000000000 

3 

0.5555555555 5555555556 
0.8888888888 8888888889 

4 

0.3478548451 3745385737 
0.6521451548 6254614263 

5 

0.2369268850 5618908751 
0.4786286704 9936646804 
0.5688888888 8888888889 

6 

0.1713244923 7917034504 
0.3607615730 4813860757 
0.4679139345 7269104739 

7 

0.1294849661 6886969327 
0.2797053914 8927666790 
0.3818300505 0511894495 
0.4179591836 7346938776 

8 

0.1012285362 9037625915 
0.2223810344 5337447054 
0.3137066458 7788728734 
0.3626837833 7836198297 


(- 1 ) 0.6667134430 8688137594 
0.1494513491 5058059315 
0.2190863625 1598204400 
0.2692667193 0999635509 
0.2955242247 1475287017 

12 

(- 1 ) 0.4717533638 6511827195 
0.1069393259 9531843096 
0.1600783285 4334622633 
0.2031674267 2306592175 
0.2334925365 3835480876 
0.2491470458 1340278500 


0.9894009349 9164993260 
0.9445750230 7323257608 
0.8656312023 8783174388 
0.7554044083 5500303390 
0.6178762444 0264374845 
0.4580167776 5722738634 
0.2816035507 7925891323 
(- 1 ) 0.9501250983 7637440185 


16 

(- 1 ) 0. 2715245941 1754094852 
(- 1 ) 0.6225352393 8647892863 
(- 1 ) 0.9515851168 2492784810 
0.1246289712 5553387205 
0.1495959888 1657673208 
0.1691565193 9500253819 
0.1826034150 4492358887 
0.1894506104 5506849629 
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